
gem = readfits('./grs_paper/DHT22_GemOB1_mom.fits',hd)
ancil = omni_read_conffile('./conffiles/ancillary.conf')
restore,'./ancillary/bessel_parallaxes.sav',/ver
restore,'./local/BGPS_V2_velocities.sav',/ver

gem[where(gem LE -0.5, nnan)] = !values.d_nan

barr = (findgen(sxpar(hd,'NAXIS3'))-sxpar(hd,'CRPIX3')+1)*sxpar(hd,'CDELT3') + $
       sxpar(hd,'CRVAL3')
bi = where(barr GE -1.0 AND barr LE 1.5)

print,n_elements(gem),long(total(~finite(gem)))

lv = transpose(total(gem[*,*,bi],3,/NAN))*abs(sxpar(hd,'CDELT3'))
bv = transpose(total(gem,2,/NAN))*abs(sxpar(hd,'CDELT2'))
lb = total(gem,1,/NAN)*abs(sxpar(hd,'CDELT1'))

slv = size(lv,/dim)
sbv = size(bv,/dim)
slb = size(lb,/dim)

sxaddpar,hd,'NAXIS',2

lvhd = hd
bvhd = hd
lbhd = hd

;; LV
sxaddpar,lvhd,'NAXIS1',slv[0]
sxaddpar,lvhd,'NAXIS2',slv[1]

sxaddpar,lvhd,'CTYPE1',sxpar(hd,'CTYPE2')
sxaddpar,lvhd,'CRVAL1',sxpar(hd,'CRVAL2')
sxaddpar,lvhd,'CDELT1',sxpar(hd,'CDELT2')
sxaddpar,lvhd,'CRPIX1',sxpar(hd,'CRPIX2')
sxaddpar,lvhd,'CTYPE2',sxpar(hd,'CTYPE1')
sxaddpar,lvhd,'CRVAL2',sxpar(hd,'CRVAL1')
sxaddpar,lvhd,'CDELT2',sxpar(hd,'CDELT1')
sxaddpar,lvhd,'CRPIX2',sxpar(hd,'CRPIX1')
sxdelpar,lvhd,'NAXIS3'
sxdelpar,lvhd,'CTYPE3'
sxdelpar,lvhd,'CRVAL3'
sxdelpar,lvhd,'CDELT3'
sxdelpar,lvhd,'CRPIX3'
sxaddpar,lvhd,'CTYPE2','GLAT-CAR'

;; LB
sxaddpar,lbhd,'NAXIS1',slb[0]
sxaddpar,lbhd,'NAXIS2',slb[1]
sxaddpar,lbhd,'CTYPE1',sxpar(hd,'CTYPE2')
sxaddpar,lbhd,'CRVAL1',sxpar(hd,'CRVAL2')
sxaddpar,lbhd,'CDELT1',sxpar(hd,'CDELT2')
sxaddpar,lbhd,'CRPIX1',sxpar(hd,'CRPIX2')
sxaddpar,lbhd,'CTYPE2',sxpar(hd,'CTYPE3')
sxaddpar,lbhd,'CRVAL2',sxpar(hd,'CRVAL3')
sxaddpar,lbhd,'CDELT2',sxpar(hd,'CDELT3')
sxaddpar,lbhd,'CRPIX2',sxpar(hd,'CRPIX3')
sxdelpar,lbhd,'NAXIS3'
sxdelpar,lbhd,'CTYPE3'
sxdelpar,lbhd,'CRVAL3'
sxdelpar,lbhd,'CDELT3'
sxdelpar,lbhd,'CRPIX3'
sxaddpar,lbhd,'CTYPE2','GLAT-CAR'


;;******************************************************************************
;; PLOT-PLOT-PLOT-PLOT-PLOT-PLOT-PLOT-PLOT-PLOT-PLOT-PLOT-PLOT
myps,'./grs_paper/gem_ob1_analyze.eps',xsize=3.1,ysize=2.8,thick=1.5
print,!x.omargin,!y.omargin
!x.omargin = [0.25,1.25]
!y.omargin = [-0.75,-2.5]
multiplot,[1,2]

csf = 0.55

undefine,pr
print,'LB PR: ',set_plot_range(lb)
pr = [2,40]                     ; K km/s
plotmap,lb,lbhd,xr=[194.1,187.9],yr=[-0.9,1.7],min_dpi=100,charsize=1.0*csf,$
        xtickformat='blank_axis',xtit='',range=pr,ct=75;,ytickformat="(F0.1)"

;; Draw in outline of BGPS coverage
cgPlots,[193.5,190.5,190.5,187.5,187.5,190.5,190.5,193.5,193.5],[0.5,0.5,1.0,1.0,0.0,0.0,-0.5,-0.5,0.5],color='wt1',linestyle=2,thick=2.5



pos = [0.90,!y.window[0],0.93,!y.window[1]]
cgColorbar,/vert,/right,range=pr,title='!u12!nCO Intensity  [K km s!u-1!n]',$
           charsize=0.9*csf,position=pos

al_legend,/top,/left,box=0,['Gem OB1'],charsize=1.0*csf

kk = where( abs(bessel.l - 191.) LE 3.0 AND abs(bessel.b) LE 1.5, nkk) 

print,'NAMES: ',1.d/bessel[kk].px

;; Right side
tenpc = 10.d / (1.d3 / bessel[kk[1]].px) / !dtor
print,'10 pc: ',tenpc
cgPlots,[188.9,188.9-tenpc],[+0.1,+0.1],thick=4
cgText,188.9-tenpc/2.,align=0.5,+0.15,'10 pc',charsize=0.7*csf

;; Left side
tenpc = 10.d / (1.d3 / bessel[kk[2]].px) / !dtor
print,'10 pc: ',tenpc
cgPlots,[193,193-tenpc],[+0.8,+0.8],thick=4
cgText,193-tenpc/2.,align=0.5,+0.85,'10 pc',charsize=0.7*csf


y0 = bessel[kk].vlsr - ancil.ppv_dv
y1 = bessel[kk].vlsr + ancil.ppv_dv

dl = ancil.ppv_dis / (1.d3 / bessel[kk].px) / !dtor
x0 = bessel[kk].l - dl
x1 = bessel[kk].l + dl

FOR pp=0,nkk-1 DO BEGIN
   tvcircle,dl[pp],bessel[kk[pp]].l,bessel[kk[pp]].b,'blk7',/data
   
   cgOplot,bessel[kk[pp]].l,bessel[kk[pp]].b,psym='filledstar',$
           color='bisque',symsize=1.2
   cgOplot,bessel[kk[pp]].l,bessel[kk[pp]].b,psym='star',$
           color='black',symsize=1.2,thick=1
ENDFOR

cgOplot,v.l,v.b,psym=16,symsize=0.2
print,'LB PR: ',pr

loff = [0.90,-0.05,-0.37]
boff = [+0.65,-0.50,+0.15]

cgText,bessel[kk].l+loff,bessel[kk].b+boff,bessel[kk].alias,$
       charsize=0.9*csf             ;,align=0.5


;; Plot N,E markers
euler,191,0.25,ra,dec,1
north = transform_posang(0,ra,dec,1)
east  = transform_posang(270,ra,dec,1)
print,'Position Angles: ',north,east

cgArrow,188.5,-0.7,188.5-0.4*sin(north*!dtor),-0.7+0.4*cos(north*!dtor),/solid,$
        color='black',/data,thick=1.5,hthick=0.25
cgArrow,188.5,-0.7,188.5-0.4*sin(east*!dtor),-0.7+0.4*cos(east*!dtor),/solid,$
        color='black',/data,thick=1.5,hthick=0.25

cgText,188.5-0.4*sin(north*!dtor)-0.05,-0.7+0.4*cos(north*!dtor)-0.05,$
       'N',charsize=0.8*csf
cgText,188.5-0.4*sin(east*!dtor)-0.1,-0.7+0.4*cos(east*!dtor),$
       'E',charsize=0.8*csf

cgAxis,yaxis=1,/yst,ytickformat='blank_axis'

multiplot,/doyaxis

;;============================================
;; L-V

undefine,pr
pr = [0,5.4] ; K arcdeg
plotmap,lv,lvhd,ct=75,axsel=7,/notsquare,yr=[-10,15],range=pr,xmargin=[10,10],$
        charsize=1.0*csf,xr=[194.1,187.9],min_dpi=100
print,'LV PR: ',pr

pos = [0.90,!y.window[0],0.93,!y.window[1]]
cgColorbar,/vert,/right,range=pr,title='!u12!nCO Intensity  [K arcdeg]',$
           charsize=0.9*csf,position=pos


FOR pp=0,nkk-1 DO BEGIN
   print,x0[pp],x1[pp],y0[pp],y1[pp]
   
   cgPolygon,[x0[pp],x1[pp],x1[pp],x0[pp],x0[pp]],$
             [y0[pp],y0[pp],y1[pp],y1[pp],y0[pp]],color='blk7',/data,thick=2.5
   cgOplot,bessel[kk[pp]].l,bessel[kk[pp]].vlsr,psym='filledstar',$
           color='bisque',symsize=1.2
   cgOplot,bessel[kk[pp]].l,bessel[kk[pp]].vlsr,psym='star',$
           color='black',symsize=1.2,thick=1
   
ENDFOR
cgOplot,v.l,v.vlsr,psym=16,symsize=0.2

loff = [+1.25,-0.65,+0.0]
voff = [-3.0,-1.0,-5.5]

cgText,bessel[kk].l+loff,bessel[kk].vlsr+voff,bessel[kk].alias,$
       charsize=0.9*csf,align=0.5,color='black'

myps,/done,/mp



END
